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Abstract — Mid-infrared astronomy (operating at wavelengths 
ranging from 2 to 25 /xm) has progressed significantly in the last 
decades, thanks to the improvement of detector techniques and 
the growing diameter of telescopes. Space observatories benefit 
from the absence of atmospheric absorption, allowing to reach 
the very high sensitivities needed to perform 3D hyperspectral 
observations, but telescopes are limited in diameter (< 1 meter) 
and therefore provide observations at low angular resolution 
(typically a few seconds of arc). On the other hand, ground- 
based facilities suffer from strong atmospheric absorption but use 
large telescopes (above 8m diameter) to perform sub-arcsecond 
angular resolution imaging through selected windows in the mid- 
infrared range. In this Paper, we present a method based on 
Lee and Seung's Non-negative Matrix Factorization (NMF) to 
merge data from space and ground based mid-infrared (mid- 
IR) telescopes in order to combine the best sensitivity, spectral 
coverage and angular resolution. We prove the efficiency of this 
technique when applied to real mid-IR astronomical data. We 
suggest that this method can be applied to any combination of 
low and high spatial resolution positive hyperspectral datasets, 
as long as the spectral variety of the data allows decomposition 
into components using NMF. 

Index Terms — Astronomy, Telescopes, Spectroscopy, Signal 
processing, Hypercubes 



I. Introduction 

In galaxies (including our own Galaxy, the Milky- Way), 
the ultraviolet (UV) and visible light (respectively wavelength 
ranges of 10-400 and 400-750 nm) emitted by stars is absorbed 
by dust particles having sizes ranging from a few nanometers 
to a few micrometers. The UV- visible energy they absorb is 
then re-emitted in the infrared (IR, 2-1000/im fTl). Therefore, 
in the recent years, astronomers have focused their efforts in 
studying the Milky-Way and external galaxies to the IR, where 
the emitted light contains information both on the amount of 
absorbed energy originating from stars and on the composition 
of interstellar dust. Unfortunately, most of the IR light coming 
from space is absorbed by the atmosphere. This has motivated 
a number of IR space missions that provide the sensitivity 
required to perform hyperspectral observations. However, due 
to technical and cost constrains, space mission can only launch 
small diameter telescopes, which limits the spatial resolution 
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of the obtained data. Conversely, the ground-based telescopes 
working in the IR usually have large apertures providing 
subarcsecond angular resolution, but low sensitivities, that 
only permit imaging through a few broad band filters. In 
the field of remote sensing, a set of methods referred to as 
pansharpening ^ have been developed in order to perform 
data fusion/improvement by combining hyperspectral datasets 
at different spatial resolutions. In this paper, we present such a 
method, developed to combine ground and space data, to ben- 
efit form the advantages of the two techniques. This method is 
based on the decomposition into components of hyperspectral 
data obtained by space telescopes, using Non-Negative Matrix 
Factorization (NMF), followed by non-negative least square 
fitting of ground-based data using these components (Sect|II|. 
We apply this to real data obtained in the mid-IR range (5-15 
/im) in order to prove the efficiency of the proposed method 
(SectHni). 

II. Formalism 

A. Decomposition into components with NMF 

We define the hyperspectral observations of a region of the 
sky as a 3 dimensional rrisXnsX Is matrix Cs(x, A) where 
(x, y) define the spatial coordinates and A the spectral index. 
We assume that all the points in Cs{x^y^ A) are positive. We 
call spectrum each vector Xs{px-,Vy, A) recorded at a position 
{PxjPy) over the Is wavelength points. We define a new 
positive 2D matrix of observations X^, the rows of which 
contain the rris x Us = ks, Xs spectra of Cs • We now assume 
that each spectrum Xs is the result of a linear combination 
of a limited number r^, with << ks, of unknown source 
spectra written ^^(A) so that: 



(1) 



where the a^{px-,Py) coefficients are unknown. This can be 
re-written under the following matrix product: 



Xs — ^.s X Ss 



(2) 



where As is the ks x matrix of unknown coefficients of the 
linear combinations and Ss is an x Is matrix, the rows of 
which are the source spectra. This is a typical blind source 
separation (BSS) problem [3J, and can be solved using multi- 
ple methods (e.g. (H, O, O). Here, we concentrate on Non- 
Negative matrix factorization lH that is applicable because A 
and S are positive. The objective is to find estimations of As 
and Ss, respectively Ws and Hs so that 



X, ^ I^, X Hs 



(3) 
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This is done by adapting the non-negative matrices Ws and 
Hs so as to minimise the squared EucHdian distance — 
Ws Hs\\^ or the divergence D{Xs\WsHs), respectively defined 
as 

- WsH4^ = yaii - {w^H^rf, (4) 



ij 



and 



(5) 
(6) 



where the exponents i and j respectively refer to the row 
and column indexes of the matrices. The algorithm used to 
achieve the minimisation of the Euclidian distance is based 
on the iterative update rule 



(7) 



and for divergence 



Ws ^ W.. 



Y,,wixi/iw,H,y 



(8) 



Euclidian distance and divergence are non increasing under 
their respective update rules, so that starting from random 
Ws and Hs matrices, the algorithm will converge towards 
a minimum for these criteria. This provides the matrix Hs 
containing the estimated source spectra. 

B. Panshapening 

Once this step has been achieved, one can observe the same 
region, with a much higher spatial resolution, but much fewer 
points in wavelength. This is typically what is done in the 
mid-IR from ground-based telescopes. This will provide a new 
hyperspectral cube ^, A) with kg = rrig x Ug » rris x 

Us spatial positions and Ig « Is spectral points. As in Sect. 
III-AI we define an observation matrix Xg that contains the 
spectra of Cg. Again we assume that each spectrum in Xg 
can be approximated using a linear combination of source 
spectra s^with Ig points. In matrix form, this read: 



Xa 



Wg X Hg, 



(9) 



where Xg is the matrix of spectra observed from the ground, 
Hg is the matrix of source spectra, and Wg a matrix of un- 
known coefficients. In fact, the source spectra over Is spectral 
points have already been estimated above using NMF and are 
given by Hs. Because the same region is being observed, the 
source spectra in Hg are expected to be the same as those 
in Hs but recorder over Ig spectral positions. Therefore, we 
construct Hg by taking the points in Hs corresponding to the 
Ig spectral positions. Now that Xg and Hg are known, Eq.[9]no 



longer defines a BSS problem and Wg can easily be adjusted 
by least square methods H. A the end of this process we 
therefore have in Wg the kg coefficients by which the Hg 
matrix has to be multiplied to recover the observed spectrum. 
Instead, one can build a new pansharped matrix Xp defined 
by: 



X^ = WaX Hs 



(10) 



This matrix contains the spectra for kg spatial positions and Is 
wavelength points. Finally, the spectra in the Xp matrix can 
be reassigned their respective positions in an rrig x rig x Is 
hyperspectral cube Cp{x^y^X). The Cp cube has therefore the 
high number of spatial points of Cg and high number of points 
in wavelength of Cg. 



III. Application to real data 

In order to test the above method we have applied it to real 
astronomical data of the NGC 7023 nebula. The observed set 
that will serve as reference (Fig. l.d) is a spectral cube Cref 
obtained by the Spitzer space Telescope |^ recorded over 36 x 
28 spatial positions and 180 points in wavelength in the mid- 
IR (5-14 /im). We simulate a low spatial resolution cube by 
degrading the spatial resolution of Cref by a factor of 4. This 
provides Cs that has 9x7 spatial points but still 180 points 
in wavelength (Fig. l.a). Cg is constructed by keeping only 
10 points in wavelengths of Cref, corresponding to windows 
that can be observed from the ground, but with the full spatial 
resolution. Cg has 36 x 28 spatial points but only 10 spectral 
positions (Fig. l.b.). The first step consists in applying NMF to 
Cs as described in Sect jII-AI In order to do this, the number 
of rows Ts of Hs has to be set. To do this we apply NMF 
for successive values of Ts and keep the smallest value that 
provides: 



{Xs-WsXHsf ^^N^ 



(11) 



where Ns represents additive contribution of noise to signal 
in Cs. We find that rg = 3 satisfies the above statement, and 
the 3 obtained spectra of Hs are presented in Fig. l.a. In the 
present application, each of the obtained source spectra can be 
attributed to the emission of different chemical species: neutral 
and ionized poly cyclic aromatic hydrocarbon molecules (resp. 
R\H^ and R\H+) and very small carbonaceous grains (VSG). 
This attribution has been discussed in [9| and ifTOl and is 
beyond the scope of the present paper. However, we will keep 
these acronyms in the following for practical reasons. Using 
the obtained Hs matrix, we construct the Hg matrix following 
the strategy explained in Sect. III-BI Using the Non Negative 
Least Square (NNLS) algorithm [7|, we identify the values of 
the coefficients in Wg. Finally, using equation [TOl we build Xp 
from Wg and Hs. We then reshape Xp in order to recover a 
spectral cube Cp with 36 x 28 spatial positions and 180 points 
(Fig. 1. c). 
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Fig. 1. Illustration of the steps achieved to perform pansharpening on mid-IR data of the NGC 7023 nebula: a) Application of non-negative matrix 
factorization to the low spatial resolution C s cube. The source spectra in Hs'. PAH*-", PAH"^ and VSG are shown in the lower panel. The colors in the upper 
image show the respective contribution, given by Ws, of these three source spectra in the C s cube at each spatial position. Colors combine as standard 
RGB i.e. green-i-red=yellow etc. b) Upper: image at high spatial resolution taken at a given wavelength position from the C g cube. Lower panel shows the 
spectra at position indicated by a cross in the image, c) High spatial/sepctral resolution C p cube obtained using NMF pansharpening. The contribution of 
each source spectrum from Hg at each spatial position is given with RGB colors as for the C s cube. Lower panel shows the reconstructed full mid-IR 
spectrum at the position indicated on the image by a cross, d) Reference cube C ref- The colors show the contribution of each source spectra from Hg to 
the mid-IR spectrum at each spatial position, obtained by NNLS. Lower panel shows the reference spectrum a the position indicated by a cross in the image. 
In all panels, the black mask in the upper right-hand corner of the images is used to mask the contribution from a bright star. 



A. Efficiency of the method 

The efficiency of the proposed method can be estimated by 
comparing Cp and Cref- Visually, the reconstruction seems 
very efficient (See Fig. 1 c)) and d)) for the spatial distribution 
of the source spectra. Fig. 2 shows an overlay of the original 
observed spectrum at the position marked with a cross in Fig. 
1 and the reconstructed spectrum. The match is excellent. 
Quantitatively, we can estimate the reconstruction efficiency 
of our method, Q, defined by: 



Q 



(12) 



which compares the norm of the residuals cube Cp — Cref 
and the norm in the reference cube Cref - We find that Q ~ 
0.09 meaning that the average reconstruction precision is 9% 
which is similar to the intrinsic uncertainty present in the Cref 



data. The full application presented above was implemented 
under Matlab. Given that this method only involves multi- 
plicative calculation, the whole process NMF+Pansharpening 
runs in less then a minute on a laptop computer. We however 
note that this method will only perform for datasets which 
show significant spectral variations, i.e. providing a variety 
of observed mixtures that is sufficient to achieve NMF. Fur- 
thermore, the number of spectral points obtained in the high 
spatial resolution dataset must be sufficient (typically about 
10) to allow a unique solution to NNLS fitting. 

B. Astrophysical implications 

We have proven the efficiency of the proposed method 
to obtain "pansharped" cubes with full spatial and spectral 
samplings, combining data from space and ground-based 
observatories. This has several advantages for astronomical 
observations: 
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Fig. 2. A 10 points spectrum taken from Cg at the position marked with 
a cross in Fig. l(red crosses). The continuous line shows the spectrum taken 
from Cp at the same position, reconstructed by fitting the ten points with a 
linear combination of source spectra. In diamonds is the reference spectrum 
taken in C^ef at the same position. 



- Astronomical spectra are always positive so this method 
can easily be applied. 

- The method is shown to increase the spatial resolution of 
mid-IR hyperspecral observations obtained by space telescopes 
(here by a factor 4). 

- The method can provide the full mid-IR spectrum, includ- 
ing in the spectral region that not accessible from the ground. 

- The proposed technique can yield spatially resolved mid- 
IR hyperspectral observations of astronomical objects for 
which it would be impossible by classical techniques. 

- The procedure involves only multiplicative computations, 
so is easy to implement and fast. 

The next step consists in obtaining real mid-IR high angular 
resolution data from a ground-based telescope, for astro- 
nomical objects previously observed with the Spitzer Space 
Telescope. We have obtained 10 hours of telescope time 
on the largest single aperture telescope in the world, the 
Grantecan (GTC) telescope, in order to do this. We will then 
be able to combine Spitzer and GTC observations using non- 
negative matrix pansharpening. This will provide the highest 
angular resolution (0.2 arceseconds) mid-IR spectral cube ever 
obtained for an astronomical object. 

IV. CONCLUSION 

We have proposed an NMF-based pansharpening method 
that allows to benefit from the high spatial resolution of mid- 
IR instruments in ground-based observatories and sensitivity 
and spectral coverage of space telescopes. The only working 
hypothesis is that the data is positive, which in the case of re- 
mote sensing is usually the case. We have successfully applied 
this technique to real astronomical observations. Promising ap- 
plications could be performed in the context of future ground- 
and space-based mid-IR missions. In particular, the NASA 
James Webb Space Telescope (JWST hereafter) and SPICA 
(JAXA) space missions will provide mid-IR spectral cubes at 
angular resolutions of 0.2 arcseconds. Meanwhile, the future 



ground-based, 42 meters, European Large Telescope (ELT 
hereafter) will observe at milliarcsecond angular resolutions. 
The combination of JWST/SPICA data with ELT data using 
NMF pansharpening will provide data at the scale size of, for 
example, the habitable zone of planet forming disks around 
young stellar objects. Finally, we emphasize the fact that NMF 
pansharpening can well be apply to any combination of low 
and high spatial resolution hyperspectral datasets for which an 
NMF decomposition into parts is found. 
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